author: “Journ Galvan” date: “4/23/2020” output: html_document —
# Load libraries
library(here)
library(dplyr)
library(tidyverse)
library(vegan)
library(reshape)
library(ggplot2)
library(ggpubr)
library(zoo)
library(sjPlot)
library(sjmisc)
library(sjlabelled)
library(psych)
library(knitr)
library(faraway)
library(car)
library(MASS)
library(gridExtra)
library(grid)
library(AER)
#Load data
branch_width <- read.csv(here("photogrammetry_data/branchwidth_data.csv"))
coral_pg <- read.csv(here("photogrammetry_data/galvan_journ_datasheet_v2.csv"))
coral_field <- read.csv(here("photogrammetry_data/field_experiment_colony_measurements_moorea_summer2019.csv"))
cafi_surveys <- read.csv(here("cafi_data/revised_cafi_data_moorea_summer2019_11_27.csv"))
cafi_field <- read.csv(here("cafi_data/prelim_cafi_counts_moorea_summer2019.csv"))
updated_cafi <- read.csv(here("cafi_data/cafi_data_w_taxonomy_summer2019_2020_5_21.csv"))
Data was filtered to only include invertebrates identified in coral colonies. These were identified down to the species level. Species richness, abundance and Shannon weiner diversity index was calculated.
#Incorporate cleaned CAFI data
updated_cafi2 <- updated_cafi %>% filter(str_detect(coral_id, "^FE")) %>%
dplyr::select(master_sort, coral_id, code, type, search_term, lowest_level, phylum, genus, species, general_notes) %>%
subset(phylum!="Chordata") #filter out vertebrates
#Calculate CAFI richness, abundance, and diversity (shannon weiner) for each coral
cafi_summarized2 <- group_by(updated_cafi2, coral_id) %>%
summarise(num_cafi = n(), cafi_richness = length(unique(code)), cafi_present = paste(sort(unique(code)), collapse = ";"))
cafi_summarized2$sw <- updated_cafi2 %>%
count(code, coral_id = coral_id) %>%
spread(code,n) %>%
mutate_all(list(~tidyr::replace_na(.,0))) %>%
dplyr::select(-coral_id) %>%
diversity(index = "shannon")
Cleaned data of coral colonies taken from the field and using photogrammetry techniques were combined into one dataframe. Many more and accurate morphological measurements could be acquired using photogrammetry than with traditional measurements taken *in situ. Common measurements of volume, surface area and available space known as interstitial space were calculated of each coral colonly.
Available space was traditionally measured using 2D measurements of five randomly selected branch distances and averaging them. In this experiment, we took methods used by Neil E. Doszpot et al. to calculate a 3D measurement of available space using convex hull geometry and the software estimated coral volume.
Because we did not take random branch distance measurements in the field, we used the software models to randomly select branch distances to be averaged.
Three other measurements were calculated from the coral models: convexity and sphericity which capture volume compactness and packing which captures how much of an objects surface area is situated internally versus externally.
Wide or tight branching corals were initially qualitatively determined. However, classifying corals with convexity values \(>= 0.5\) as “tight” and \(< 0.5\) as “wide” offers another method of classifying coral morphology. For the rest of these tests, we will classify “wide” and “tight” branching corals using this method.
# Clean field data
coral_field2 <- coral_field %>%
dplyr::rename(branch = branch_width)
# Clean photogrammetry morphometric data
coral_pg$volume_pg <- as.numeric(as.character(coral_pg$volume_pg)) # Convert factor to numeric
# Clean photogrammetry branch width data and take averages of branch distance for each coral
branch_width$branch_distance_mm <- as.numeric(as.character(branch_width$branch_distance_mm)) # Convert factor to numeric
branch_w_summarized <- group_by(branch_width, coral_id) %>%
summarise(avg_w_cm = sum(branch_distance_mm/10)/n(), #take average branch distance and convert mm to cm
measurements = length(unique(replicate_measurement)),
locations = paste(sort(unique(location)), collapse = ";"))
#Create dataframe called coral_dim
coral_dim <- merge(coral_pg, branch_w_summarized, by = "coral_id") %>%
merge(coral_field2, by = "coral_id") %>%
mutate(volume_pg=volume_pg*10^6,
max_hull_volume=max_hull_volume*10^6, #convert m^3 to cm^3
max_hull_surface_area=max_hull_surface_area*10^4,
surface_area=surface_area*10^4,
height_pg=height_pg*100,
length_pg=length_pg*100,
width_pg=width_pg*100) %>% #convert m^2 to cm^2
dplyr::rename(max_hull_SA = max_hull_surface_area,
SA = surface_area)
coral_dim$interstitial_space <- coral_dim$max_hull_volume - coral_dim$volume_pg #calculate available space by subtracting software estimated volume from convex hull volume
coral_dim$SAV <- coral_dim$SA / coral_dim$volume_pg #calculate surface area to volume relationship
coral_dim$convexity <- coral_dim$volume_pg / coral_dim$max_hull_volume #calculate proportion occupied, high ratios indicate less free space between branches
coral_dim$packing <- coral_dim$SA / coral_dim$max_hull_SA #how much of an objects surface area is situated internally
coral_dim$sphericity <- ((3.14^(1/3))*((6*coral_dim$volume_pg)^(2/3)))/coral_dim$SA #calculate sphericity or how close the object is to a sphere
coral_dim$diff_est <- coral_dim$volume_field-coral_dim$volume_pg #actual overestimation
coral_dim$prop_est <- (coral_dim$volume_field - coral_dim$volume_pg)/coral_dim$volume_pg #proportion overestimated
coral_dim <- coral_dim %>%
dplyr::select(coral_id, branch, cafi, size_class, volume_field, volume_pg, SA, height_pg, length_pg, width_pg, avg_w_cm, interstitial_space, SAV, convexity, packing, sphericity, diff_est, prop_est) %>%
filter(cafi=="empty")
#coral_dim$branch <- ifelse(coral_dim$convexity>=0.5, "tight","wide") #classifies wide and tight branching coral based on convexity measurement
Dataframes were combined and seperate columns created for log and square root transformations.
#Create dataframe called cafi_coral and merge coral and cafi data
cafi_coral <- merge(cafi_summarized2, coral_dim, by ="coral_id") %>%
drop_na(interstitial_space) #%>%
#dplyr::select(coral_id, branch, cafi, size_class, volume_field, volume_pg, SA, height_pg, length_pg, width_pg, avg_w_cm, interstitial_space, SAV, convexity, packing, sphericity, diff_est, prop_est, num_cafi, cafi_present, cafi_richness, sw)
cafi_coral$branch <- ifelse(cafi_coral$convexity>=0.5, "tight","wide") #classifies wide and tight branching coral based on convexity measurement
#Log transform variables
cafi_coral$log_num_cafi <- log(cafi_coral$num_cafi)
cafi_coral$log_volume_field <- log(cafi_coral$volume_field)
cafi_coral$log_volume_pg <- log(cafi_coral$volume_pg)
cafi_coral$log_avg_w <- log(cafi_coral$avg_w_cm)
cafi_coral$log_interstitial_space <- log(cafi_coral$interstitial_space)
cafi_coral$log_SA <- log(cafi_coral$SA)
cafi_coral$log_cafi_richness <- log(cafi_coral$cafi_richness)
cafi_coral$log_sphericity <- log(cafi_coral$sphericity)
#Sqrt transform variables
cafi_coral$sqrt_num_cafi <- sqrt(cafi_coral$num_cafi)
cafi_coral$sqrt_volume_pg <- sqrt(cafi_coral$volume_pg)
cafi_coral$sqrt_avg_w_cm <- sqrt(cafi_coral$avg_w_cm)
cafi_coral$sqrt_volume_field <- sqrt(cafi_coral$volume_field)
cafi_coral$sqrt_interstitial_space <- sqrt(cafi_coral$interstitial_space)
cafi_coral$sqrt_SA <- sqrt(cafi_coral$SA)
cafi_coral$sqrt_cafi_richness <- sqrt(cafi_coral$cafi_richness)
cafi_coral$sqrt_sphericity <- sqrt(cafi_coral$sphericity)
Look for unimodal or bimodal curve in data for average branch distance, interstitial space and convexity. Compare these curves to branch distances defined in the field to branch distances defined by convexity measuremnt of 0.5
#Load plyr package
library(plyr)
## ------------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## ------------------------------------------------------------------------------
##
## Attaching package: 'plyr'
## The following object is masked from 'package:faraway':
##
## ozone
## The following object is masked from 'package:ggpubr':
##
## mutate
## The following objects are masked from 'package:reshape':
##
## rename, round_any
## The following object is masked from 'package:purrr':
##
## compact
## The following objects are masked from 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
## The following object is masked from 'package:here':
##
## here
#Classification based on convexity measuremnts
mu <- ddply(cafi_coral, "branch",
summarise, grp.mean=mean(avg_w_cm, na.rm=TRUE))
d1 <- ggplot(cafi_coral, aes(x=avg_w_cm, fill=branch))+
geom_density(alpha=0.4)+
geom_vline(data=mu, aes(xintercept=grp.mean, color=branch),
linetype="dashed")+
labs(x="Branch distance")+
theme(legend.position = "top")
mu2 <- ddply(cafi_coral, "branch",
summarise, grp.mean=mean(interstitial_space, na.rm=TRUE))
d2 <- ggplot(cafi_coral, aes(x=interstitial_space, fill=branch))+
geom_density(alpha=0.4)+
geom_vline(data=mu2, aes(xintercept=grp.mean, color=branch),
linetype="dashed")+
labs(x="Interstitial space")+
theme(legend.position = "top")
mu3 <- ddply(cafi_coral, "branch",
summarise, grp.mean=mean(convexity, na.rm=TRUE))
d3 <- ggplot(cafi_coral, aes(x=convexity, fill=branch))+
labs(x="Convexity")+
geom_density(alpha=0.4)+
geom_vline(data=mu3, aes(xintercept=grp.mean, color=branch),
linetype="dashed")+
theme(legend.position = "top")
#Classification based on qualitative analysis in the field
mu4 <- ddply(coral_dim, "branch",
summarise, grp.mean=mean(avg_w_cm, na.rm=TRUE))
d4 <- ggplot(coral_dim, aes(x=avg_w_cm, fill=branch))+
geom_density(alpha=0.4)+
geom_vline(data=mu4, aes(xintercept=grp.mean, color=branch),
linetype="dashed")+
labs(x="Branch distance")+
theme(legend.position = "top")
mu5 <- ddply(coral_dim, "branch",
summarise, grp.mean=mean(interstitial_space, na.rm=TRUE))
d5 <- ggplot(coral_dim, aes(x=interstitial_space, fill=branch))+
geom_density(alpha=0.4)+
geom_vline(data=mu5, aes(xintercept=grp.mean, color=branch),
linetype="dashed")+
labs(x="Interstitial space")+
theme(legend.position = "top")
mu6 <- ddply(coral_dim, "branch",
summarise, grp.mean=mean(convexity, na.rm=TRUE))
d6 <- ggplot(coral_dim, aes(x=convexity, fill=branch))+
geom_density(alpha=0.4)+
geom_vline(data=mu6, aes(xintercept=grp.mean, color=branch),
linetype="dashed")+
labs(x="Convexity")+
theme(legend.position = "top")
fig1 <- ggarrange(d1,d2,d3,d4,d5,d6,
ncol = 3, nrow = 2)
annotate_figure(fig1,
fig.lab = "Figure 1", fig.lab.face = "bold")
Non-transformed and transformed morphological measurements correlated using *pairs.panels().
#Visualize correlations between variables
cor_var <- cafi_coral %>%
dplyr::select(num_cafi,volume_pg,height_pg,volume_field,SA,avg_w_cm,interstitial_space,sphericity,convexity,packing)
pairs.panels(cor_var,
method="pearson",
hist.col = "white",
density=TRUE,
ellipses=FALSE,
scale=TRUE)
#Use sqrt transformed data
cor_var2 <- cafi_coral %>%
dplyr::select(sqrt_num_cafi,sqrt_volume_pg,sqrt_volume_field,sqrt_SA,sqrt_avg_w_cm,sqrt_interstitial_space,sqrt_sphericity,convexity,packing)
pairs.panels(cor_var2,
method="pearson",
hist.col = "white",
density=TRUE,
ellipses=FALSE,
scale=TRUE)
#Use log transformed data
cor_var3 <- cafi_coral %>%
dplyr::select(log_num_cafi,log_volume_pg,log_volume_field,log_SA,log_avg_w,log_interstitial_space,log_sphericity,convexity,packing)
pairs.panels(cor_var3,
method="pearson",
hist.col = "white",
density=TRUE,
ellipses=FALSE,
scale=TRUE)
Log transformed field and software estimated volume and log transformed software estimated SA and interstitial space are all highly correlated to eachother. Packing is also highly correlated to volume and SA.
Here, we want to compare the differences between our software and manual measurements of volume, height, length and width. Manual measurements for coral volume were estimated from an elipsoid.
dim_pg <- coral_pg %>%
dplyr::select(coral_id, height_pg, length_pg, width_pg, volume_pg) %>%
mutate(height_pg = height_pg*100,
length_pg = length_pg*100,
width_pg = width_pg*100,
volume_pg = volume_pg*10^6)
dim_field <- coral_field %>%
dplyr::select(coral_id, height_field, length_field, width_field, volume_field)
dim_bind <- merge(dim_pg, dim_field, by="coral_id") %>%
drop_na() %>%
summarise(coral_id,volume_pg,volume_field,
height_pg,height_field,
length_pg,length_field,
width_pg,width_field,
prop_est=(volume_field-volume_pg)/volume_pg, #proportion overestimation of elipsoide compared to software volume
diff_est=volume_field-volume_pg) #actual overestimation difference
Since the corals are not independent, that is, they were taken on the same coral but with different methods, we will consider a paired t-test with dependent samples. Height, length and width measurements are also subjective since they may have been taken from different areas on the coral. We will first check for normallity assumptions. Normality will be checked with a qqplot.
par(mfrow=c(3,3))
qqPlot(dim_bind$height_pg, dist="norm")
## [1] 1 58
qqPlot(dim_bind$height_field, dist="norm")
## [1] 1 56
qqPlot(dim_bind$length_pg, dist="norm")
## [1] 6 57
qqPlot(dim_bind$length_field, dist="norm")
## [1] 6 57
qqPlot(dim_bind$width_pg, dist="norm")
## [1] 19 1
qqPlot(dim_bind$width_field, dist="norm")
## [1] 1 23
par(mfrow=c(2,2))
par(mfrow=c(2,2))
qqPlot(dim_bind$volume_pg, dist="norm")
## [1] 58 23
qqPlot(dim_bind$volume_field, dist="norm")
## [1] 1 36
par(mfrow=c(1,1))
#Volume data is not normal, try sqrt transformation
dim_bind$sqrt_volume_pg <- sqrt(dim_bind$volume_pg)
dim_bind$sqrt_volume_field <- sqrt(dim_bind$volume_field)
par(mfrow=c(2,2))
qqPlot(dim_bind$sqrt_volume_pg, dist="norm")
## [1] 58 23
qqPlot(dim_bind$sqrt_volume_field, dist="norm")
## [1] 1 36
par(mfrow=c(1,1))
#sqrt transformed volume for both measurements is normal, move to paired t-test
t.test(dim_bind$sqrt_volume_pg,dim_bind$sqrt_volume_field, paired=TRUE)
##
## Paired t-test
##
## data: dim_bind$sqrt_volume_pg and dim_bind$sqrt_volume_field
## t = -8.8976, df = 58, p-value = 1.957e-12
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -29.08742 -18.40334
## sample estimates:
## mean of the differences
## -23.74538
#ddply(dim_bind, .(sqrt_volume_pg,sqrt_volume_field), summarize,
t.test(dim_bind$height_pg,dim_bind$height_field, paired=TRUE)
##
## Paired t-test
##
## data: dim_bind$height_pg and dim_bind$height_field
## t = -5.6455, df = 58, p-value = 5.194e-07
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -2.483686 -1.183436
## sample estimates:
## mean of the differences
## -1.833561
t.test(dim_bind$length_pg,dim_bind$length_field, paired=TRUE)
##
## Paired t-test
##
## data: dim_bind$length_pg and dim_bind$length_field
## t = -0.0097688, df = 58, p-value = 0.9922
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.5147722 0.5097722
## sample estimates:
## mean of the differences
## -0.0025
t.test(dim_bind$width_pg,dim_bind$width_field, paired=TRUE)
##
## Paired t-test
##
## data: dim_bind$width_pg and dim_bind$width_field
## t = -4.2154, df = 58, p-value = 8.836e-05
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -1.4710549 -0.5237823
## sample estimates:
## mean of the differences
## -0.9974186
We reject the null hypothesis that the true difference in means is equal to zero for all measurements except for length.
#Log SA
g1 <- ggplot(cafi_coral, aes(x=log_volume_pg, y=log_SA, shape=branch, color=branch, linetype=branch))+
geom_point()+
geom_smooth(method=lm, se=TRUE, aes(fill=branch))+
theme_classic()+
theme(axis.title.x=element_blank(),
legend.position="none")+
labs(y=expression(paste("Log SA cm"^{2})))
#Log avg. width
g2 <- ggplot(cafi_coral, aes(x=log_volume_pg, y=log_avg_w, shape=branch, color=branch, linetype=branch))+
geom_point()+
geom_smooth(method=lm, se=TRUE, aes(fill=branch))+
theme_classic()+
theme(axis.title.x=element_blank(),
legend.position="none")+
labs(y="Log branch distance cm")
#Log IS
g3 <- ggplot(cafi_coral, aes(x=log_volume_pg, y=log_interstitial_space, shape=branch, color=branch, linetype=branch))+
geom_point()+
geom_smooth(method=lm, se=TRUE, aes(fill=branch))+
theme_classic()+
theme(axis.title.x=element_blank(),
legend.position="none")+
labs(y=expression(paste("Log space cm"^{3})))
#Log Sphericity
g4 <- ggplot(cafi_coral, aes(x=log_volume_pg, y=log_sphericity, shape=branch, color=branch, linetype=branch))+
geom_point()+
geom_smooth(method=lm, se=TRUE, aes(fill=branch))+
theme_classic()+
theme(axis.title.x=element_blank(),
legend.position="none")+
labs(y="Log sphericity")
#Convexity
g5 <- ggplot(cafi_coral, aes(x=log_volume_pg, y=convexity, shape=branch, color=branch, linetype=branch))+
geom_point()+
geom_smooth(method=lm, se=TRUE, aes(fill=branch))+
theme_classic()+
theme(axis.title.x=element_blank(),
legend.position="none")+
labs(y="Convexity", x=expression(paste("Log Vol. cm"^{3})))
#Packing
g6 <- ggplot(cafi_coral, aes(x=log_volume_pg, y=packing, shape=branch, color=branch, linetype=branch))+
geom_point()+
geom_smooth(method=lm, se=TRUE, aes(fill=branch))+
theme_classic()+
theme(axis.title.x=element_blank(),
legend.position="bottom")+
labs(y="Packing")
#extract legend
#https://github.com/hadley/ggplot2/wiki/Share-a-legend-between-two-ggplot2-graphs
g_legend<-function(a.gplot){
tmp <- ggplot_gtable(ggplot_build(a.gplot))
leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
legend <- tmp$grobs[[leg]]
return(legend)}
mylegend<-g_legend(g6)
fig2 <- ggarrange(g1,g2,g3,g4,g5,g6,
ncol = 3, nrow=2)
annotate_figure(fig2,
top = text_grob("Morphological Measurements", color = "black", face = "bold", size = 14),
bottom = text_grob(expression(paste("Log Vol. cm"^{3})), color = "black",
face = "bold", size = 12),
fig.lab = "Figure 2", fig.lab.face = "bold")
The below graphs compare the software estimated volume, height, length and width measurements to the same coral’s manual measurements. A one to one line was fit to show how differences in measuremets vary with increasing coral size.
#Volume
#summary(lm(formula = volume_field ~ volume_pg, data = dim_bind))
volume <- ggplot(dim_bind, aes(x=volume_pg, y=volume_field))+
geom_point()+
geom_smooth(method=lm, se=TRUE)+
scale_y_continuous(breaks=seq(0,max(dim_bind$volume_field), by = 5000))+
scale_x_continuous(breaks=seq(0,max(dim_bind$volume_pg), by = 2000))+
theme_bw()+
labs(y= expression(paste("elipsoid V cm"^{3})), x= expression(paste("software V cm"^{3})))+
geom_abline(intercept = 756.3176, slope = 1, color='black') #x-intercept pulled from summary data
#Height
#summary(lm(formula = height_field ~ height_pg, data = dim_bind))
height <- ggplot(dim_bind, aes(x=height_pg, y=height_field))+
geom_point()+
geom_smooth(method=lm, se=TRUE)+
theme_bw()+
labs(y= "Manual h cm", x= "Software h cm")+
geom_abline(intercept = 2.44284, slope = 1, color='black')
#Length
#summary(lm(formula = length_field ~ length_pg, data = dim_bind))
length <- ggplot(dim_bind, aes(x=length_pg, y=length_field))+
geom_point()+
geom_smooth(method=lm, se=TRUE)+
theme_bw()+
labs(y= "Manual l cm", x= "Software l cm")+
geom_abline(intercept = 0.52596, slope = 1, color='black')
#Width
#summary(lm(formula = width_field ~ width_pg, data = dim_bind))
width <- ggplot(dim_bind, aes(x=width_pg, y=width_field))+
geom_point()+
geom_smooth(method=lm, se=TRUE)+
theme_bw()+
labs(y= "Manual w cm", x= "Software w cm")+
geom_abline(intercept = 2.13398, slope = 1, color='black')
fig3 <- ggarrange(volume, height, length, width,
ncol = 2, nrow = 2)
annotate_figure(fig3,
top = text_grob("Manual v. software measurements", color = "black", face = "bold", size = 14),
fig.lab = "Figure 3", fig.lab.face = "bold")
The below graphs show the proportion and actual differences of software volume from that of elipsoid volume. This change is shown with increasing software coral volume which we will assume is the closest to the actual volume of the coral.
#Actual overestimation
plot_est <- ggplot(cafi_coral, aes(x=volume_pg, y=diff_est, shape=branch, color=branch, linetype=branch))+
geom_point()+
geom_smooth(method=lm, se=TRUE, aes(fill=branch))+
theme_bw()+
labs(y= expression(paste("elipsoid overestim. cm"^{3})), x= expression(paste("software V cm"^{3})))+
theme(legend.position="none")
#Linear Model
lm_plot_est <- lm(diff_est~volume_pg+branch, data=cafi_coral)
#Check normality, homoscedasticity, linearity
par(mfrow=c(2,2))
plot(lm_plot_est)
par(mfrow=c(1,1))
#Proportional overestimation
plot_prop <- ggplot(cafi_coral, aes(x=volume_pg, y=prop_est, shape=branch, color=branch, linetype=branch))+
geom_point()+
geom_smooth(method=lm, se=TRUE, aes(fill=branch))+
theme_bw()+
labs(y= "Proportion overestim.", x= expression(paste("software V cm"^{3})))
#Linear Model
lm_plot_prop <- lm(prop_est~volume_pg+branch, data=cafi_coral)
#Check normality, homoscedasticity, linearity
par(mfrow=c(2,2))
plot(lm_plot_prop)
par(mfrow=c(1,1))
fig1 <- ggarrange(plot_est, plot_prop,
ncol = 2, nrow = 1)
annotate_figure(fig1,
top = text_grob("Actual vs. proportional estimates", color = "black", face = "bold", size = 14),
fig.lab = "Figure 1", fig.lab.face = "bold")
ggsave("fig1.pdf", plot = fig1, path="figs", width=7, height=4)
#Create table with the linear model parameters
pl <- c(
'(Intercept)' = "Intercept",
volume_pg = "Software V")
pl2 <- c(
'(Intercept)' = "Intercept",
volume_pg = "Software V")
tab_model(lm_plot_est,lm_plot_prop,
pred.labels = c(pl,pl2),
show.se = TRUE,
string.ci = "Conf. Int (95%)",
string.p = "P-Value",
dv.labels = c("Overestim.", "Proportion over-est."))
| Overestim. | Proportion over-est. | |||||||
|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | std. Error | Conf. Int (95%) | P-Value | Estimates | std. Error | Conf. Int (95%) | P-Value |
| Intercept | -2587.16 | 2075.50 | -6853.42 – 1679.09 | 0.224 | 0.51 | 0.27 | -0.05 – 1.07 | 0.075 |
| Software V | 1.14 | 0.26 | 0.61 – 1.67 | <0.001 | 0.00 | 0.00 | -0.00 – 0.00 | 0.333 |
| branchwide | 4843.95 | 2066.56 | 596.07 – 9091.82 | 0.027 | 1.06 | 0.27 | 0.50 – 1.62 | 0.001 |
| Observations | 29 | 29 | ||||||
| R2 / R2 adjusted | 0.434 / 0.390 | 0.371 / 0.323 | ||||||
Measuring the volume of an elipsoid overestimates actual coral volume approximately by 50% across all size ranges and overestimates this volume by a difference of approximately \(4280.195 cm^3\) Looking at the plots, as our best estimate of coral volume increases, the difference between the elipsoid volume and our best estimate also increases in a positive relationship. From our second plot of overestimation proportion with increasing coral volume as the predictor variable, there is no evidence of the proportion or \(software V/elipsoid V\) increasing with increasing best estimate volume with a \(p-value=0.061\). The proportion stayed constant at approximately 50% overestimation of our best estimate of the actual coral volume. ########################################################################################################################
The below graphs compare the sqrt and log transformed software estimated volume to the same coral’s sqrt and log transformed elipsoid volume. A one to one line was fit to show how differences in measuremets vary with increasing coral size.
#Sqrt volume
#summary(lm(formula = sqrt_volume_field ~ sqrt_volume_pg, data = cafi_coral))
sqrt_volume <- ggplot(cafi_coral, aes(x=sqrt_volume_pg, y=sqrt_volume_field))+
geom_point()+
geom_smooth(method=lm, se=TRUE)+
theme_bw()+
labs(y= expression(paste("sqrt elipsoid V cm"^{3})), x= expression(paste("sqrt software V cm"^{3})))+
geom_abline(intercept = 4.207, slope = 1, color='blue') #x-intercept pulled from summary data
#summary(lm(formula = log_volume_field ~ log_volume_pg, data = cafi_coral))
log_volume <- ggplot(cafi_coral, aes(x=log_volume_pg, y=log_volume_field))+
geom_point()+
geom_smooth(method=lm, se=TRUE)+
theme_bw()+
labs(y= expression(paste("log elipsoid V cm"^{3})), x= expression(paste("log software V cm"^{3})))+
geom_abline(intercept = 0.80673, slope = 1, color='blue') #x-intercept pulled from summary data
fig4 <- ggarrange(volume, sqrt_volume, log_volume,
ncol = 3, nrow = 1)
annotate_figure(fig4,
top = text_grob("Compare transformations", color = "black", face = "bold", size = 14),
fig.lab = "Figure 4", fig.lab.face = "bold")
The below graphs represent how invertebrate abundance and richness relates to both volume measurements. We will test how well our software estimated method of calculating volume and interstitial space relates to invertebrate counts. We will first assess \(p<0.05\) to test if the linear regression slope was significantly different from zero. Correlation of dependent and independent variables \(R^2\) was also included in the tables.
#Plot invert abundance to elipsoid volume
invert_field <- ggplot(cafi_coral, aes(x=volume_field, y=num_cafi))+
geom_point()+
geom_smooth(method=lm, se=TRUE)+
theme_classic()+
labs(y="invert abundance", x= expression(paste("elipsoid V cm"^{3})))
#Check qqplots for normality
lm_invert_field <- lm(num_cafi~volume_field, data=cafi_coral)
par(mfrow=c(2,2))
plot(lm_invert_field)
par(mfrow=c(1,1))
#Plot invert abundance to software estimated volume
invert_pg <- ggplot(cafi_coral, aes(x=volume_pg, y=num_cafi))+
geom_point()+
geom_smooth(method=lm, se=TRUE)+
theme_classic()+
labs(y="invert abundance", x= expression(paste("software V cm"^{3})))
#Check qqplots for normality
lm_invert_pg <- lm(num_cafi~volume_pg, data=cafi_coral)
par(mfrow=c(2,2))
plot(lm_invert_pg)
par(mfrow=c(1,1))
#Plot invert richness to software estimated volume
richness_field <- ggplot(cafi_coral, aes(x=volume_field, y=cafi_richness))+
geom_point()+
geom_smooth(method=lm, se=TRUE)+
theme_classic()+
labs(y="invert richness", x= expression(paste("elipsoid V cm"^{3})))
#Check qqplots for normality
lm_richness_field <- lm(cafi_richness~volume_field, data=cafi_coral)
par(mfrow=c(2,2))
plot(lm_richness_field)
par(mfrow=c(1,1))
#Plot invert richness to software estimated volume
richness_pg <- ggplot(cafi_coral, aes(x=volume_pg, y=cafi_richness))+
geom_point()+
geom_smooth(method=lm, se=TRUE)+
theme_classic()+
labs(y="invert richness", x= expression(paste("software V cm"^{3})))
#Check qqplots for normality
lm_richness_pg <- lm(cafi_richness~volume_pg, data=cafi_coral)
par(mfrow=c(2,2))
plot(lm_richness_pg)
par(mfrow=c(1,1))
fig2 <- ggarrange(invert_field, invert_pg, richness_field, richness_pg,
ncol = 2, nrow = 2)
annotate_figure(fig2,
top = text_grob("Invert Abundance and Richness", color = "black", face = "bold", size = 14),
fig.lab = "Figure 2", fig.lab.face = "bold")
ggsave("fig2.pdf", plot = fig2, path="figs", width=5, height=4)
#HTML table comparing volume methods and abundance
pl3 <- c(
'(Intercept)' = "Intercept",
volume_field = "Elipsoid volume")
pl4 <- c(
'(Intercept)' = "Intercept",
volume_pg = "Software est. volume")
tab_model(lm_invert_field,lm_invert_pg,
pred.labels = c(pl3,pl4),
show.se = TRUE,
string.ci = "Conf. Int (95%)",
string.p = "P-Value",
digits=5,
dv.labels = c("Field volume vs. abundance","Software volume vs. abundance"))
| Field volume vs. abundance | Software volume vs. abundance | |||||||
|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | std. Error | Conf. Int (95%) | P-Value | Estimates | std. Error | Conf. Int (95%) | P-Value |
| Intercept | 13.69118 | 3.08747 | 7.35622 – 20.02615 | <0.001 | 11.54653 | 2.92781 | 5.53916 – 17.55390 | 0.001 |
| Elipsoid volume | 0.00052 | 0.00024 | 0.00002 – 0.00102 | 0.043 | ||||
| Software est. volume | 0.00169 | 0.00053 | 0.00060 – 0.00277 | 0.004 | ||||
| Observations | 29 | 29 | ||||||
| R2 / R2 adjusted | 0.143 / 0.111 | 0.275 / 0.248 | ||||||
#HTML table comparing volume methods and richness
pl5 <- c(
'(Intercept)' = "Intercept",
volume_field = "Elipsoid volume")
pl6 <- c(
'(Intercept)' = "Intercept",
volume_pg = "Software est. volume")
tab_model(lm_richness_field,lm_richness_pg,
pred.labels = c(pl5,pl6),
show.se = TRUE,
string.ci = "Conf. Int (95%)",
string.p = "P-Value",
digits=5,
dv.labels = c("Field volume vs. richness","Software volume vs. richness"))
| Field volume vs. richness | Software volume vs. richness | |||||||
|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | std. Error | Conf. Int (95%) | P-Value | Estimates | std. Error | Conf. Int (95%) | P-Value |
| Intercept | 6.24931 | 0.85828 | 4.48826 – 8.01037 | <0.001 | 5.87648 | 0.84250 | 4.14782 – 7.60514 | <0.001 |
| Elipsoid volume | 0.00019 | 0.00007 | 0.00005 – 0.00033 | 0.011 | ||||
| Software est. volume | 0.00051 | 0.00015 | 0.00019 – 0.00082 | 0.003 | ||||
| Observations | 29 | 29 | ||||||
| R2 / R2 adjusted | 0.218 / 0.189 | 0.291 / 0.265 | ||||||
When comparing elipsoid volume to software estimated volume for invertebrate richness, we see significant p-values for both as represented in the above tables. Software estimated volume measurements and invert. richness show a lower p-value indicating a greater slope from zero as well as having a higher correlation \(R^2\) value.
Software estimated available space calculated from convex hull and coral volume was related to invertebrate abundance. The same was repeated for the traditional method of average branch distance.
Dependent and independent variables were also compared after being log transformed.
branch_pg <- ggplot(cafi_coral, aes(x=interstitial_space, y=num_cafi))+
geom_point()+
geom_smooth(method=lm, se=TRUE)+
theme_bw()+
labs(y="invert abundance", x= expression(paste("IS cm"^{3})))
lm_branch_pg <- lm(num_cafi~interstitial_space, data=cafi_coral)
log_branch_pg <- ggplot(cafi_coral, aes(x=log_interstitial_space, y=log_num_cafi))+
geom_point()+
geom_smooth(method=lm, se=TRUE)+
theme_bw()+
labs(y="Log invert abundance", x= expression(paste("Log IS cm"^{3})))
lm_log_branch_pg <- lm(log_num_cafi~log_interstitial_space, data=cafi_coral)
branch_field <- ggplot(cafi_coral, aes(x=avg_w_cm, y=num_cafi))+
geom_point()+
geom_smooth(method=lm, se=TRUE)+
theme_bw()+
labs(y="invert abundance", x= "avg_w_cm")
lm_branch_field <- lm(num_cafi~avg_w_cm, data=cafi_coral)
log_branch_field <- ggplot(cafi_coral, aes(x=log_avg_w, y=log_num_cafi))+
geom_point()+
geom_smooth(method=lm, se=TRUE)+
theme_bw()+
labs(y="Log invert abundance", x= "Log avg. width")
lm_log_branch_field <- lm(log_num_cafi~log_avg_w, data=cafi_coral)
fig7 <- ggarrange(branch_field, branch_pg, log_branch_field, log_branch_pg,
ncol = 2, nrow = 2)
annotate_figure(fig7,
top = text_grob("Invert abundance and habitat space", color = "black", face = "bold", size = 14),
fig.lab = "Figure 7", fig.lab.face = "bold")
#HTML table from regression results
pl <- c(
'(Intercept)' = "Intercept",
interstitial_space = "IS")
tab_model(lm_branch_pg,
pred.labels = pl,
dv.labels = "Invertebrate abundance")
| Invertebrate abundance | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| Intercept | 16.01 | 9.85 – 22.17 | <0.001 |
| IS | 0.00 | -0.00 – 0.00 | 0.250 |
| Observations | 29 | ||
| R2 / R2 adjusted | 0.049 / 0.013 | ||
#HTML table from regression results
pl <- c(
'(Intercept)' = "Intercept",
log_interstitial_space = "Log IS")
tab_model(lm_log_branch_pg,
pred.labels = pl,
dv.labels = "Invertebrate abundance")
| Invertebrate abundance | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| Intercept | 0.56 | -0.59 – 1.70 | 0.329 |
| Log IS | 0.28 | 0.13 – 0.43 | 0.001 |
| Observations | 29 | ||
| R2 / R2 adjusted | 0.354 / 0.330 | ||
#HTML table from regression results
pl <- c(
'(Intercept)' = "Intercept",
avg_w_cm = "Avg. branch width")
tab_model(lm_branch_field,
pred.labels = pl,
dv.labels = "Invertebrate abundance")
| Invertebrate abundance | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| Intercept | 17.06 | 5.45 – 28.66 | 0.006 |
| Avg. branch width | 0.60 | -5.02 – 6.21 | 0.829 |
| Observations | 29 | ||
| R2 / R2 adjusted | 0.002 / -0.035 | ||
#HTML table from regression results
pl <- c(
'(Intercept)' = "Intercept",
log_avg_w = "Log avg. width")
tab_model(lm_log_branch_field,
pred.labels = pl,
dv.labels = "Invertebrate abundance")
| Invertebrate abundance | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| Intercept | 2.54 | 2.10 – 2.98 | <0.001 |
| Log avg. width | 0.23 | -0.42 – 0.87 | 0.476 |
| Observations | 29 | ||
| R2 / R2 adjusted | 0.019 / -0.017 | ||
Log transformed data appeared to not violate normality assumptions so was used in this analysis. Relating interstitial space to invertebrate abundance gave a signficant p-value while average width did not. The difference in p-values and \(R^2\) correlation values were significant. This may provide evidence to average branch width not adequately capturing available space for invertebrates when compared to a 3D software estimated measurement.
branch_rich_pg <- ggplot(cafi_coral, aes(x=interstitial_space, y=cafi_richness))+
geom_point()+
geom_smooth(method=lm, se=TRUE)+
theme_bw()+
labs(y="invert richness", x= expression(paste("IS cm"^{3})))
lm_branch_rich_pg <- lm(cafi_richness~interstitial_space, data=cafi_coral)
log_branch_rich_pg <- ggplot(cafi_coral, aes(x=log_interstitial_space, y=log_cafi_richness))+
geom_point()+
geom_smooth(method=lm, se=TRUE)+
theme_bw()+
labs(y="Log invert richness", x= expression(paste("Log IS cm"^{3})))
lm_log_branch_rich_pg <- lm(log_cafi_richness~log_interstitial_space, data=cafi_coral)
branch_rich_field <- ggplot(cafi_coral, aes(x=avg_w_cm, y=cafi_richness))+
geom_point()+
geom_smooth(method=lm, se=TRUE)+
theme_bw()+
labs(y="invert richness", x= "avg_w_cm")
lm_branch_rich_field <- lm(cafi_richness~avg_w_cm, data=cafi_coral)
log_branch_rich_field <- ggplot(cafi_coral, aes(x=log_avg_w, y=log_cafi_richness))+
geom_point()+
geom_smooth(method=lm, se=TRUE)+
theme_bw()+
labs(y="Log invert richness", x= "Log avg. width")
lm_log_branch_rich_field <- lm(log_cafi_richness~log_avg_w, data=cafi_coral)
fig8 <- ggarrange(branch_rich_field, branch_rich_pg, log_branch_rich_field, log_branch_rich_pg,
ncol = 2, nrow = 2)
annotate_figure(fig8,
top = text_grob("Invert richness and habitat space", color = "black", face = "bold", size = 14),
fig.lab = "Figure 8", fig.lab.face = "bold")
#HTML table from regression results
pl <- c(
'(Intercept)' = "Intercept",
interstitial_space = "IS")
tab_model(lm_branch_rich_pg,
pred.labels = pl,
dv.labels = "Invertebrate richness")
| Invertebrate richness | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| Intercept | 6.83 | 5.11 – 8.54 | <0.001 |
| IS | 0.00 | -0.00 – 0.00 | 0.053 |
| Observations | 29 | ||
| R2 / R2 adjusted | 0.131 / 0.099 | ||
#HTML table from regression results
pl <- c(
'(Intercept)' = "Intercept",
log_interstitial_space = "Log IS")
tab_model(lm_log_branch_rich_pg,
pred.labels = pl,
dv.labels = "Invertebrate richness")
| Invertebrate richness | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| Intercept | 0.33 | -0.47 – 1.14 | 0.403 |
| Log IS | 0.21 | 0.11 – 0.32 | <0.001 |
| Observations | 29 | ||
| R2 / R2 adjusted | 0.390 / 0.368 | ||
#HTML table from regression results
pl <- c(
'(Intercept)' = "Intercept",
avg_w_cm = "Avg. branch width")
tab_model(lm_branch_rich_field,
pred.labels = pl,
dv.labels = "Invertebrate richness")
| Invertebrate richness | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| Intercept | 6.50 | 3.17 – 9.83 | <0.001 |
| Avg. branch width | 0.73 | -0.88 – 2.34 | 0.361 |
| Observations | 29 | ||
| R2 / R2 adjusted | 0.031 / -0.005 | ||
#HTML table from regression results
pl <- c(
'(Intercept)' = "Intercept",
log_avg_w = "Log avg. width")
tab_model(lm_log_branch_rich_field,
pred.labels = pl,
dv.labels = "Invertebrate richness")
| Invertebrate richness | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| Intercept | 1.81 | 1.50 – 2.13 | <0.001 |
| Log avg. width | 0.25 | -0.22 – 0.71 | 0.285 |
| Observations | 29 | ||
| R2 / R2 adjusted | 0.042 / 0.007 | ||
Log transformed data appeared to not violate normality assumptions so was used in this analysis. Relating interstitial space to invertebrate richness gave a signficant p-value while average width did not. The difference in p-values and \(R^2\) correlation values were significant. This may provide evidence to average branch width not adequately capturing available space for invertebrates when compared to a 3D software estimated measurement.
Visualize data using boxplots.
# Boxplot of invert abundance
log_invert_box <- ggplot(cafi_coral, aes(x= branch, y=log_num_cafi, color = branch))+
geom_boxplot(width = 0.5)+
geom_jitter()+ #adds data on top of box plot
scale_color_manual(values = c("firebrick1", "steelblue2"))+
scale_y_continuous(limits= c(0, max(cafi_coral$log_num_cafi + 1)))+
theme_classic()+
labs(y="Log invert. abundance", x="Coral branching", title = "Log invert. abundance boxplot")
log_invert_box
# Boxplot of invertebrate species richness per coral for both tight and wide branching
richness_box <- ggplot(cafi_coral, aes(x= branch, y=log_cafi_richness, color = branch))+
geom_boxplot(width = 0.5)+
geom_jitter()+
scale_color_manual(values = c("firebrick1", "steelblue2"))+
scale_y_continuous(limits= c(0, max(cafi_coral$log_cafi_richness + 1)))+
theme_classic()+
labs(y="Invert. sp. richness", x="Coral branching", title = "Invert. sp. richness boxplot")
richness_box
# Boxplot of invertebrate diversity for both tight and wide branching coral
diversity_box <- ggplot(cafi_coral, aes(x= branch, y=sw, color = branch))+
geom_boxplot(width = 0.5)+
geom_jitter()+
scale_color_manual(values = c("firebrick1", "steelblue2"))+
scale_y_continuous(limits= c(0, max(cafi_coral$sw + 0.5)), breaks = c(0, 0.5, 1, 1.5, 2, 2.5, 3))+
theme_classic()+
labs(y="Invert. diversity", x="Coral branching", title = "Invert. diversity boxplot")
diversity_box
Data was log transformed except for diversity so as not to violate normality assumptions. It appears that there is no difference in abundance, species richness or diversity.
Before conducting a two sample t-test we checked for normality using qq-plots and Shapiro Wilk test and a Levene’s test for equal variance.
#Q-Q plot separated by group
with(cafi_coral, qqPlot(num_cafi[branch == "wide"], dist="norm"))
## [1] 13 14
with(cafi_coral, qqPlot(num_cafi[branch == "tight"], dist="norm"))
## [1] 12 6
#Shapiro Wilk seperated by group
with(cafi_coral, shapiro.test(num_cafi[branch == "wide"]))
##
## Shapiro-Wilk normality test
##
## data: num_cafi[branch == "wide"]
## W = 0.92775, p-value = 0.1995
with(cafi_coral, shapiro.test(num_cafi[branch == "tight"]))
##
## Shapiro-Wilk normality test
##
## data: num_cafi[branch == "tight"]
## W = 0.76166, p-value = 0.003547
With a \(p>0.05\) we fail to reject the null hypothesis that our data is normally distributed for wide branching corals. With a \(p<0.05\) we reject the null hypothesis that our data is normally distributed for tight branching corals. Because we are working with smaller data sets, we cannot use the central limit theorem.
#Levene's test: response variable, group variable
leveneTest(cafi_coral$num_cafi, cafi_coral$branch)
## Warning in leveneTest.default(cafi_coral$num_cafi, cafi_coral$branch):
## cafi_coral$branch coerced to factor.
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.2017 0.6569
## 27
With a \(p>0.05\) we fail to reject the null hypothesis that our data is normally distributed for tight and wide branching corals. Both groups of corals that have tight and wide branching morphologies have equal variance.
Because tight branching coral is not normally distributed, we could either perform a transformation or non-parametric test. We will try log transforming our data first. A column log transforming the number of inverts has already been created
#Q-Q plot separated by group
with(cafi_coral, qqPlot(log_num_cafi[branch == "wide"], dist="norm"))
## [1] 8 6
with(cafi_coral, qqPlot(log_num_cafi[branch == "tight"], dist="norm"))
## [1] 12 6
#Shapiro Wilk seperated by group
with(cafi_coral, shapiro.test(log_num_cafi[branch == "wide"]))
##
## Shapiro-Wilk normality test
##
## data: log_num_cafi[branch == "wide"]
## W = 0.93506, p-value = 0.2642
with(cafi_coral, shapiro.test(log_num_cafi[branch == "tight"]))
##
## Shapiro-Wilk normality test
##
## data: log_num_cafi[branch == "tight"]
## W = 0.88255, p-value = 0.09449
#Levene's test: response variable, group variable
leveneTest(cafi_coral$log_num_cafi, cafi_coral$branch)
## Warning in leveneTest.default(cafi_coral$log_num_cafi, cafi_coral$branch):
## cafi_coral$branch coerced to factor.
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.4491 0.5084
## 27
# subset data
tight <- subset(cafi_coral, branch == "tight")
wide <- subset(cafi_coral, branch == "wide")
t.test(tight$log_num_cafi, wide$log_num_cafi, var.equal = TRUE)
##
## Two Sample t-test
##
## data: tight$log_num_cafi and wide$log_num_cafi
## t = 0.9573, df = 27, p-value = 0.3469
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.3005325 0.8262331
## sample estimates:
## mean of x mean of y
## 2.81695 2.55410
With a \(p>0.05\) we fail to reject the null hypothesis that the difference between the means of the log transformed invert counts are different from zero.
####Compare other photogrametry morphological measurements to invert abundance
#Volume
plot1 <- ggplot(cafi_coral, aes(x=volume_pg, y=num_cafi))+
geom_point()+
geom_smooth(method=lm, se=TRUE)+
theme_classic()+
theme(axis.title.y=element_blank(),
legend.position="none")+
labs(x=expression(paste("Volume cm"^{3})))
#Sphericity
plot2 <- ggplot(cafi_coral, aes(x=sphericity, y=num_cafi))+
geom_point()+
geom_smooth(method=lm, se=TRUE)+
theme_classic()+
theme(axis.title.y=element_blank(),
legend.position="none")+
labs(x="Sphericity")
lm_sphericity <- lm(num_cafi~sphericity, data = cafi_coral)
#Packing
plot3 <- ggplot(cafi_coral, aes(x=packing, y=num_cafi))+
geom_point()+
geom_smooth(method=lm, se=TRUE)+
theme_classic()+
theme(axis.title.y=element_blank(),
legend.position="none")+
labs(x="Packing")
lm_packing <- lm(num_cafi~packing, data = cafi_coral)
#Height
plot4 <- ggplot(cafi_coral, aes(x=height_pg, y=num_cafi))+
geom_point()+
geom_smooth(method=lm, se=TRUE)+
theme_classic()+
theme(axis.title.y=element_blank(),
legend.position="none")+
labs(x="Height")
lm_height <- lm(num_cafi~height_pg, data = cafi_coral)
#Length
plot5 <- ggplot(cafi_coral, aes(x=length_pg, y=num_cafi))+
geom_point()+
geom_smooth(method=lm, se=TRUE)+
theme_classic()+
theme(axis.title.y=element_blank(),
legend.position="none")+
labs(x="Length")
lm_length <- lm(num_cafi~length_pg, data = cafi_coral)
#Width
plot6 <- ggplot(cafi_coral, aes(x=width_pg, y=num_cafi))+
geom_point()+
geom_smooth(method=lm, se=TRUE)+
theme_classic()+
theme(axis.title.y=element_blank(),
legend.position="none")+
labs(x="Width")
lm_width <- lm(num_cafi~width_pg, data = cafi_coral)
#Figure 3
fig3 <- ggarrange(plot1,plot2,plot3,plot4,plot5,plot6,
ncol = 3, nrow=2)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
annotate_figure(fig3,
top = text_grob("Morphological Measurements", color = "black", face = "bold", size = 14),
left = text_grob("Invert. Abundance", color = "black",
rot = 90, face = "bold", size = 12),
fig.lab = "Figure 3", fig.lab.face = "bold")
ggsave("fig3.pdf", plot = fig3, path="figs", width=6, height=4)
#Data tables for volume, sphericity, packing
pl7 <- c('(Intercept)' = "Intercept",
volume_pg = "Software Volume")
pl8 <- c('(Intercept)' = "Intercept",
sphericity = "Sphericity")
pl9 <- c('(Intercept)' = "Intercept",
packing = "Packing")
tab_model(lm_invert_pg,lm_sphericity, lm_packing,
pred.labels = c(pl7,pl8,pl9),
show.se = TRUE,
string.ci = "Conf. Int (95%)",
string.p = "P-Value",
digits=5,
dv.labels = c("Abundance vs. Software Volume","Abundance vs. Sphericity","Abundance vs. Packing"))
| Abundance vs. Software Volume | Abundance vs. Sphericity | Abundance vs. Packing | ||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | std. Error | Conf. Int (95%) | P-Value | Estimates | std. Error | Conf. Int (95%) | P-Value | Estimates | std. Error | Conf. Int (95%) | P-Value |
| Intercept | 11.54653 | 2.92781 | 5.53916 – 17.55390 | 0.001 | 33.93199 | 9.99889 | 13.41597 – 54.44801 | 0.002 | -18.68154 | 12.72981 | -44.80096 – 7.43788 | 0.154 |
| Software Volume | 0.00169 | 0.00053 | 0.00060 – 0.00277 | 0.004 | ||||||||
| Sphericity | -51.45528 | 31.75640 | -116.61402 – 13.70347 | 0.117 | ||||||||
| Packing | 19.19597 | 6.53838 | 5.78032 – 32.61162 | 0.007 | ||||||||
| Observations | 29 | 29 | 29 | |||||||||
| R2 / R2 adjusted | 0.275 / 0.248 | 0.089 / 0.055 | 0.242 / 0.214 | |||||||||
#Data tables for height, length, width
pl10 <- c('(Intercept)' = "Intercept",
height_pg = "Height")
pl11 <- c('(Intercept)' = "Intercept",
length_pg = "Length")
pl12 <- c('(Intercept)' = "Intercept",
width_pg = "Width")
tab_model(lm_height,lm_length, lm_width,
pred.labels = c(pl10,pl11,pl12),
show.se = TRUE,
string.ci = "Conf. Int (95%)",
string.p = "P-Value",
digits=5,
dv.labels = c("Abundance vs. Height","Abundance vs. Length","Abundance vs. Width"))
| Abundance vs. Height | Abundance vs. Length | Abundance vs. Width | ||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | std. Error | Conf. Int (95%) | P-Value | Estimates | std. Error | Conf. Int (95%) | P-Value | Estimates | std. Error | Conf. Int (95%) | P-Value |
| Intercept | 2.04637 | 4.54287 | -7.27482 – 11.36757 | 0.656 | 3.52788 | 5.19889 | -7.13936 – 14.19512 | 0.503 | 3.93804 | 5.02293 | -6.36815 – 14.24423 | 0.440 |
| Height | 0.97632 | 0.24875 | 0.46593 – 1.48671 | 0.001 | ||||||||
| Length | 0.56273 | 0.18292 | 0.18742 – 0.93804 | 0.005 | ||||||||
| Width | 0.66226 | 0.21264 | 0.22596 – 1.09856 | 0.004 | ||||||||
| Observations | 29 | 29 | 29 | |||||||||
| R2 / R2 adjusted | 0.363 / 0.340 | 0.260 / 0.232 | 0.264 / 0.237 | |||||||||
ggsave("fig3.pdf", plot = fig3, path="figs", width=6, height=4)
####Compare other photogrametry morphological measurements to invert richness
#Volume
plot7 <- ggplot(cafi_coral, aes(x=volume_pg, y=cafi_richness))+
geom_point()+
geom_smooth(method=lm, se=TRUE)+
theme_classic()+
theme(axis.title.y=element_blank(),
legend.position="none")+
labs(x=expression(paste("Volume cm"^{3})))
#Sphericity
plot8 <- ggplot(cafi_coral, aes(x=sphericity, y=cafi_richness))+
geom_point()+
geom_smooth(method=lm, se=TRUE)+
theme_classic()+
theme(axis.title.y=element_blank(),
legend.position="none")+
labs(x="Sphericity")
lm_sph_rich <- lm(cafi_richness~sphericity, data=cafi_coral)
#Packing
plot9 <- ggplot(cafi_coral, aes(x=packing, y=cafi_richness))+
geom_point()+
geom_smooth(method=lm, se=TRUE)+
theme_classic()+
theme(axis.title.y=element_blank(),
legend.position="none")+
labs(x="Packing")
lm_pack_rich <- lm(cafi_richness~packing, data=cafi_coral)
#Height
plot10 <- ggplot(cafi_coral, aes(x=height_pg, y=cafi_richness))+
geom_point()+
geom_smooth(method=lm, se=TRUE)+
theme_classic()+
theme(axis.title.y=element_blank(),
legend.position="none")+
labs(x="Height")
lm_height_rich <- lm(cafi_richness~height_pg, data=cafi_coral)
#Length
plot11 <- ggplot(cafi_coral, aes(x=length_pg, y=cafi_richness))+
geom_point()+
geom_smooth(method=lm, se=TRUE)+
theme_classic()+
theme(axis.title.y=element_blank(),
legend.position="none")+
labs(x="Length")
lm_length_rich <- lm(cafi_richness~length_pg, data=cafi_coral)
#Width
plot12 <- ggplot(cafi_coral, aes(x=width_pg, y=cafi_richness))+
geom_point()+
geom_smooth(method=lm, se=TRUE)+
theme_classic()+
theme(axis.title.y=element_blank(),
legend.position="none")+
labs(x="Width")
lm_width_rich <- lm(cafi_richness~width_pg, data=cafi_coral)
#Figure 4
fig4 <- ggarrange(plot7,plot8,plot9,plot10,plot11,plot12,
ncol = 3, nrow=2)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
annotate_figure(fig4,
top = text_grob("Morphological Measurements", color = "black", face = "bold", size = 14),
left = text_grob("Invert. Richness", color = "black",
rot = 90, face = "bold", size = 12),
fig.lab = "Figure 4", fig.lab.face = "bold")
ggsave("fig4.pdf", plot = fig4, path="figs", width=6, height=4)
#Data tables for volume, sphericity, packing
pl7 <- c('(Intercept)' = "Intercept",
volume_pg = "Software Volume")
pl8 <- c('(Intercept)' = "Intercept",
sphericity = "Sphericity")
pl9 <- c('(Intercept)' = "Intercept",
packing = "Packing")
tab_model(lm_richness_pg,lm_sphericity, lm_packing,
pred.labels = c(pl7,pl8,pl9),
show.se = TRUE,
string.ci = "Conf. Int (95%)",
string.p = "P-Value",
digits=5,
dv.labels = c("Richness vs. Software Volume","Richness vs. Sphericity","Richness vs. Packing"))
| Richness vs. Software Volume | Richness vs. Sphericity | Richness vs. Packing | ||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | std. Error | Conf. Int (95%) | P-Value | Estimates | std. Error | Conf. Int (95%) | P-Value | Estimates | std. Error | Conf. Int (95%) | P-Value |
| Intercept | 5.87648 | 0.84250 | 4.14782 – 7.60514 | <0.001 | 33.93199 | 9.99889 | 13.41597 – 54.44801 | 0.002 | -18.68154 | 12.72981 | -44.80096 – 7.43788 | 0.154 |
| Software Volume | 0.00051 | 0.00015 | 0.00019 – 0.00082 | 0.003 | ||||||||
| Sphericity | -51.45528 | 31.75640 | -116.61402 – 13.70347 | 0.117 | ||||||||
| Packing | 19.19597 | 6.53838 | 5.78032 – 32.61162 | 0.007 | ||||||||
| Observations | 29 | 29 | 29 | |||||||||
| R2 / R2 adjusted | 0.291 / 0.265 | 0.089 / 0.055 | 0.242 / 0.214 | |||||||||
#Data tables for height, length, width
pl10 <- c('(Intercept)' = "Intercept",
height_pg = "Height")
pl11 <- c('(Intercept)' = "Intercept",
length_pg = "Length")
pl12 <- c('(Intercept)' = "Intercept",
width_pg = "Width")
tab_model(lm_height_rich,lm_length_rich, lm_width_rich,
pred.labels = c(pl10,pl11,pl12),
show.se = TRUE,
string.ci = "Conf. Int (95%)",
string.p = "P-Value",
digits=5,
dv.labels = c("Richness vs. Height","Richness vs. Length","Richness vs. Width"))
| Richness vs. Height | Richness vs. Length | Richness vs. Width | ||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | std. Error | Conf. Int (95%) | P-Value | Estimates | std. Error | Conf. Int (95%) | P-Value | Estimates | std. Error | Conf. Int (95%) | P-Value |
| Intercept | 2.70318 | 1.24142 | 0.15600 – 5.25037 | 0.038 | 3.05376 | 1.43923 | 0.10070 – 6.00681 | 0.043 | 3.25370 | 1.39842 | 0.38439 – 6.12301 | 0.028 |
| Height | 0.31233 | 0.06797 | 0.17286 – 0.45181 | <0.001 | ||||||||
| Length | 0.18476 | 0.05064 | 0.08087 – 0.28866 | 0.001 | ||||||||
| Width | 0.21441 | 0.05920 | 0.09294 – 0.33587 | 0.001 | ||||||||
| Observations | 29 | 29 | 29 | |||||||||
| R2 / R2 adjusted | 0.439 / 0.418 | 0.330 / 0.305 | 0.327 / 0.302 | |||||||||
We will now fit a poisson general linear model to our invertebrate count data
poiss_field <- glm(num_cafi~volume_field, family=poisson, data=cafi_coral) #Poisson regression for invert count data using elipsoide volume
summary(poiss_field)
##
## Call:
## glm(formula = num_cafi ~ volume_field, family = poisson, data = cafi_coral)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.2927 -1.8563 -0.4155 1.3963 5.3379
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.671e+00 6.145e-02 43.465 < 2e-16 ***
## volume_field 2.340e-05 3.912e-06 5.982 2.21e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 229.18 on 28 degrees of freedom
## Residual deviance: 197.50 on 27 degrees of freedom
## AIC: 332.46
##
## Number of Fisher Scoring iterations: 5
poiss_pg <- glm(num_cafi~volume_pg, family=poisson, data=cafi_coral) #Poisson regression for invert count data using software est. volume
summary(poiss_pg)
##
## Call:
## glm(formula = num_cafi ~ volume_pg, family = poisson, data = cafi_coral)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -5.5429 -2.1822 -0.2023 1.2582 5.1645
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.544e+00 6.653e-02 38.233 <2e-16 ***
## volume_pg 7.738e-05 9.381e-06 8.248 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 229.18 on 28 degrees of freedom
## Residual deviance: 167.64 on 27 degrees of freedom
## AIC: 302.61
##
## Number of Fisher Scoring iterations: 5
The residual deviance shows overdispersion for both models. Consequently, we will run a negative binomial model to account for overdispersion. To test if elispoid and software estimated models are different, we will perform a Vuong test for non-nested models.
glm.nb1 <- glm.nb(num_cafi~volume_field, data=cafi_coral) #negative binomial model for elipsoid volume and abundance to account for overdispersion
summary(glm.nb1)
##
## Call:
## glm.nb(formula = num_cafi ~ volume_field, data = cafi_coral,
## init.theta = 3.248837444, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.88272 -0.76793 0.00909 0.44732 1.80535
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.548e+00 1.555e-01 16.383 < 2e-16 ***
## volume_field 3.594e-05 1.185e-05 3.033 0.00242 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(3.2488) family taken to be 1)
##
## Null deviance: 36.215 on 28 degrees of freedom
## Residual deviance: 29.666 on 27 degrees of freedom
## AIC: 217.94
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 3.249
## Std. Err.: 0.980
##
## 2 x log-likelihood: -211.938
glm.nb2 <- glm.nb(num_cafi~volume_pg, data=cafi_coral) #negative binomial model for software estimated volume and abundance to account for overdispersion
summary(glm.nb2)
##
## Call:
## glm.nb(formula = num_cafi ~ volume_pg, data = cafi_coral, init.theta = 3.741397829,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.18550 -0.97664 -0.05352 0.47083 1.88660
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.489e+00 1.525e-01 16.325 < 2e-16 ***
## volume_pg 8.935e-05 2.635e-05 3.391 0.000696 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(3.7414) family taken to be 1)
##
## Null deviance: 40.591 on 28 degrees of freedom
## Residual deviance: 29.749 on 27 degrees of freedom
## AIC: 214.64
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 3.74
## Std. Err.: 1.18
##
## 2 x log-likelihood: -208.637
#Data table of results
pl_nb <- c(
'(Intercept)' = "Intercept",
volume_field = "Elipsoid volume")
pl2_nb <- c(
'(Intercept)' = "Intercept",
volume_pg = "Software est. volume")
tab_model(glm.nb1,glm.nb2,
pred.labels = c(pl_nb,pl2_nb),
show.se = TRUE,
show.aicc = TRUE,
string.ci = "Conf. Int (95%)",
string.p = "P-Value",
digits=5,
dv.labels = c("Negative binomial abundance"))
| Negative binomial abundance | ||||||||
|---|---|---|---|---|---|---|---|---|
| Predictors | Incidence Rate Ratios | std. Error | Conf. Int (95%) | P-Value | Incidence Rate Ratios | std. Error | Conf. Int (95%) | P-Value |
| Intercept | 12.78321 | 0.15554 | 9.12736 – 18.08719 | <0.001 | 12.05323 | 0.15249 | 8.85505 – 16.55706 | <0.001 |
| Elipsoid volume | 1.00004 | 0.00001 | 1.00001 – 1.00007 | 0.002 | ||||
| Software est. volume | 1.00009 | 0.00003 | 1.00003 – 1.00015 | 0.001 | ||||
| Observations | 29 | 29 | ||||||
| R2 Nagelkerke | 0.283 | 0.414 | ||||||
| AIC | 218.898 | 215.597 | ||||||
#Vuongtest for abundance negative binomial models
library(nonnest2)
## This is nonnest2 0.5-5.
## nonnest2 has not been tested with all combinations of model classes.
vuongtest(glm.nb1,glm.nb2) #compare negative binomial models
##
## Model 1
## Class: negbin
## Call: glm.nb(formula = num_cafi ~ volume_field, data = cafi_coral, ...
##
## Model 2
## Class: negbin
## Call: glm.nb(formula = num_cafi ~ volume_pg, data = cafi_coral, init.theta = 3.741397829, ...
##
## Variance test
## H0: Model 1 and Model 2 are indistinguishable
## H1: Model 1 and Model 2 are distinguishable
## w2 = 0.138, p = 0.108
##
## Non-nested likelihood ratio test
## H0: Model fits are equal for the focal population
## H1A: Model 1 fits better than Model 2
## z = -0.827, p = 0.796
## H1B: Model 2 fits better than Model 1
## z = -0.827, p = 0.2042
We will now fit a poisson general linear model to our invertebrate richness data
poiss_field_r <- glm(cafi_richness~volume_field, family=poisson, data=cafi_coral) #Poisson regression for richness using elipsoide volume
summary(poiss_field_r)
##
## Call:
## glm(formula = cafi_richness ~ volume_field, family = poisson,
## data = cafi_coral)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.14226 -0.96518 -0.04043 0.72953 2.43666
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.869e+00 9.294e-02 20.116 < 2e-16 ***
## volume_field 2.004e-05 6.110e-06 3.279 0.00104 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 49.269 on 28 degrees of freedom
## Residual deviance: 39.591 on 27 degrees of freedom
## AIC: 154
##
## Number of Fisher Scoring iterations: 4
poiss_pg_r <- glm(cafi_richness~volume_pg, family=poisson, data=cafi_coral) #Poisson regression for invert richness using software est. volume
summary(poiss_pg_r)
##
## Call:
## glm(formula = cafi_richness ~ volume_pg, family = poisson, data = cafi_coral)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2600 -0.9414 0.1766 0.5776 2.4931
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.815e+00 9.864e-02 18.396 < 2e-16 ***
## volume_pg 5.619e-05 1.486e-05 3.783 0.000155 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 49.269 on 28 degrees of freedom
## Residual deviance: 36.051 on 27 degrees of freedom
## AIC: 150.46
##
## Number of Fisher Scoring iterations: 4
The residual deviance shows overdispersion for both models. Consequently, we will run a negative binomial model to account for overdispersion. To test if elispoid and software estimated models are different, we will perform a Vuong test for non-nested models.
glm.nb3 <- glm.nb(cafi_richness~volume_field, data=cafi_coral) #negative binomial model for elipsoid volume and richness to account for overdispersion
summary(glm.nb3)
##
## Call:
## glm.nb(formula = cafi_richness ~ volume_field, data = cafi_coral,
## init.theta = 20.91987063, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.91079 -0.83442 -0.01783 0.60680 2.03014
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.855e+00 1.090e-01 17.025 < 2e-16 ***
## volume_field 2.152e-05 7.547e-06 2.852 0.00435 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(20.9199) family taken to be 1)
##
## Null deviance: 35.986 on 28 degrees of freedom
## Residual deviance: 28.700 on 27 degrees of freedom
## AIC: 154.19
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 20.9
## Std. Err.: 19.5
##
## 2 x log-likelihood: -148.194
glm.nb4 <- glm.nb(cafi_richness~volume_pg, data=cafi_coral) #negative binomial model for software estimated volume and richness to account for overdispersion
summary(glm.nb4)
##
## Call:
## glm.nb(formula = cafi_richness ~ volume_pg, data = cafi_coral,
## init.theta = 33.32422167, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9972 -0.8639 0.1459 0.5140 2.2036
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.810e+00 1.090e-01 16.602 < 2e-16 ***
## volume_pg 5.725e-05 1.703e-05 3.361 0.000776 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(33.3242) family taken to be 1)
##
## Null deviance: 39.936 on 28 degrees of freedom
## Residual deviance: 29.257 on 27 degrees of freedom
## AIC: 151.72
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 33.3
## Std. Err.: 45.3
##
## 2 x log-likelihood: -145.723
#Data table of results
pl3_nb <- c(
'(Intercept)' = "Intercept",
volume_field = "Elipsoid volume")
pl4_nb <- c(
'(Intercept)' = "Intercept",
volume_pg = "Software est. volume")
tab_model(glm.nb3,glm.nb4,
pred.labels = c(pl3_nb,pl4_nb),
show.se = TRUE,
show.aicc = TRUE,
string.ci = "Conf. Int (95%)",
string.p = "P-Value",
digits=5,
dv.labels = c("Negative binomial richness"))
| Negative binomial richness | ||||||||
|---|---|---|---|---|---|---|---|---|
| Predictors | Incidence Rate Ratios | std. Error | Conf. Int (95%) | P-Value | Incidence Rate Ratios | std. Error | Conf. Int (95%) | P-Value |
| Intercept | 6.39297 | 0.10897 | 5.13025 – 7.92216 | <0.001 | 6.10986 | 0.10902 | 4.91180 – 7.54485 | <0.001 |
| Elipsoid volume | 1.00002 | 0.00001 | 1.00001 – 1.00004 | 0.004 | ||||
| Software est. volume | 1.00006 | 0.00002 | 1.00002 – 1.00009 | 0.001 | ||||
| Observations | 29 | 29 | ||||||
| R2 Nagelkerke | 0.313 | 0.412 | ||||||
| AIC | 155.154 | 152.683 | ||||||
#Vuong test for richness negative binomial models
vuongtest(glm.nb3,glm.nb4)
##
## Model 1
## Class: negbin
## Call: glm.nb(formula = cafi_richness ~ volume_field, data = cafi_coral, ...
##
## Model 2
## Class: negbin
## Call: glm.nb(formula = cafi_richness ~ volume_pg, data = cafi_coral, ...
##
## Variance test
## H0: Model 1 and Model 2 are indistinguishable
## H1: Model 1 and Model 2 are distinguishable
## w2 = 0.117, p = 0.233
##
## Non-nested likelihood ratio test
## H0: Model fits are equal for the focal population
## H1A: Model 1 fits better than Model 2
## z = -0.670, p = 0.748
## H1B: Model 2 fits better than Model 1
## z = -0.670, p = 0.2515
#Plot GLM for abundance
plot_nb1 <- ggplot(cafi_coral, aes(x=volume_field, y=num_cafi))+
geom_point()+
theme_classic()+
geom_smooth(method="glm.nb", se=TRUE)+
labs(y="Invert. Abundance", x= expression(paste("Elipsoid Volume cm"^{3})))
plot_nb2 <- ggplot(cafi_coral, aes(x=volume_pg, y=num_cafi))+
geom_point()+
theme_classic()+
geom_smooth(method="glm.nb", se=TRUE)+
labs(y="Invert. Abundance", x= expression(paste("Software Volume cm"^{3})))
#Plot GLM for richness
plot_nb3 <- ggplot(cafi_coral, aes(x=volume_field, y=cafi_richness))+
geom_point()+
theme_classic()+
geom_smooth(method="glm.nb", se=TRUE)+
labs(y="Invert. Richness", x= expression(paste("Elipsoid Volume cm"^{3})))
plot_nb4 <- ggplot(cafi_coral, aes(x=volume_pg, y=cafi_richness))+
geom_point()+
theme_classic()+
geom_smooth(method="glm.nb", se=TRUE)+
labs(y="Invert. Richness", x= expression(paste("Software Volume cm"^{3})))
#Create fig. 5 and annotate
fig5 <- ggarrange(plot_nb1, plot_nb2, plot_nb3, plot_nb4,
ncol = 2, nrow = 2)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
annotate_figure(fig5,
top = text_grob("Negative binomial", color = "black", face = "bold", size = 14),
fig.lab = "Figure 5", fig.lab.face = "bold")
ggsave("fig5.pdf", plot = fig5, path="figs", width=5, height=4)
#Create data tables from negative binomial results for abundance
pl_nb <- c(
'(Intercept)' = "Intercept",
volume_field = "Elipsoid Volume")
pl2_nb <- c(
'(Intercept)' = "Intercept",
volume_pg = "Software Volume")
tab_model(glm.nb1,glm.nb2,
pred.labels = c(pl_nb,pl2_nb),
show.se = TRUE,
show.aicc = TRUE,
string.ci = "Conf. Int (95%)",
string.p = "P-Value",
digits=5,
dv.labels = c("Negative Binomial Abundance"))
| Negative Binomial Abundance | ||||||||
|---|---|---|---|---|---|---|---|---|
| Predictors | Incidence Rate Ratios | std. Error | Conf. Int (95%) | P-Value | Incidence Rate Ratios | std. Error | Conf. Int (95%) | P-Value |
| Intercept | 12.78321 | 0.15554 | 9.12736 – 18.08719 | <0.001 | 12.05323 | 0.15249 | 8.85505 – 16.55706 | <0.001 |
| Elipsoid Volume | 1.00004 | 0.00001 | 1.00001 – 1.00007 | 0.002 | ||||
| Software Volume | 1.00009 | 0.00003 | 1.00003 – 1.00015 | 0.001 | ||||
| Observations | 29 | 29 | ||||||
| R2 Nagelkerke | 0.283 | 0.414 | ||||||
| AIC | 218.898 | 215.597 | ||||||
#Create data tables from negative binomial results for richness
pl3_nb <- c(
'(Intercept)' = "Intercept",
volume_field = "Elipsoid Volume")
pl4_nb <- c(
'(Intercept)' = "Intercept",
volume_pg = "Software Volume")
tab_model(glm.nb3,glm.nb4,
pred.labels = c(pl3_nb,pl4_nb),
show.se = TRUE,
show.aicc = TRUE,
string.ci = "Conf. Int (95%)",
string.p = "P-Value",
digits=5,
dv.labels = c("Negative Binomial Richness"))
| Negative Binomial Richness | ||||||||
|---|---|---|---|---|---|---|---|---|
| Predictors | Incidence Rate Ratios | std. Error | Conf. Int (95%) | P-Value | Incidence Rate Ratios | std. Error | Conf. Int (95%) | P-Value |
| Intercept | 6.39297 | 0.10897 | 5.13025 – 7.92216 | <0.001 | 6.10986 | 0.10902 | 4.91180 – 7.54485 | <0.001 |
| Elipsoid Volume | 1.00002 | 0.00001 | 1.00001 – 1.00004 | 0.004 | ||||
| Software Volume | 1.00006 | 0.00002 | 1.00002 – 1.00009 | 0.001 | ||||
| Observations | 29 | 29 | ||||||
| R2 Nagelkerke | 0.313 | 0.412 | ||||||
| AIC | 155.154 | 152.683 | ||||||
#Probability that the "ith" model minimizes info loss
library(AICcmodavg)
## Warning: no function found corresponding to methods exports from 'raster' for:
## 'wkt'
aicc_field <- AICc(glm.nb1) #calculate AICc from field and software models
aicc_pg <- AICc(glm.nb2)
#AIC_min
p=exp((aicc_pg-aicc_field)/2)
p
## [1] 0.1919075
With \(p=0.108\) and \(p=0.233\), we fail to reject the null hypothesis that the two models comparing elipsoid and software volume are indistinguishable for both abundance and richness using the Vuong test comparing non-nested models.
Theoretically, the lower AIC values for software estimated volume estimating log-likelihood as shown in the tables reflect a relatively better quality model to fit the data. Because a likelihood ratio test requires a nested model assumption, we used the equation \(p=exp^{AIC_{min} - AIC_i/2}\) to calculate the relative likelihood of the two models. Based on our calculation, the invertebrate model with elipsoid volume as the predictor variable was 0.19 times as probable as the model with software volume to minimize information loss.